Spatial Extremes
Whitney Huang
May 24, 2022
Acknowledgment: This R session is adopted from Ribatet’s short course on “Modelling spatial extremes with the SpatialExtremes package”.
Spatial GEV
Plot the spatial locations
library(SpatialExtremes)
data(wind)
par(mar = rep(0, 4))
maps::map(xlim = c(0, 9), ylim = c(47.5, 56))
points(coord, pch = 15, cex = 0.5)
Fitting spatial GEV
loc.form <- y ~ lon * lat; scale.form <- shape.form <- y ~ 1
#Convert coordinates from degrees to km (rough)
coord[, 1:2] <- 111 * coord[, 1:2]
(M0 <- fitspatgev(wind, scale(coord, scale=FALSE), loc.form, scale.form, shape.form))## Model: Spatial GEV model
## Deviance: 10733.34
## TIC: 10766.49
##
## Location Parameters:
## locCoeff1 locCoeff2 locCoeff3 locCoeff4
## 2.680e+02 -1.466e-01 1.888e-01 4.548e-05
## Scale Parameters:
## scaleCoeff1
## 33.99
## Shape Parameters:
## shapeCoeff1
## -0.09664
##
## Standard Errors
## locCoeff1 locCoeff2 locCoeff3 locCoeff4 scaleCoeff1 shapeCoeff1
## 3.1549723 0.0094921 0.0159158 0.0001288 1.0377212 0.0181140
##
## Asymptotic Variance Covariance
## locCoeff1 locCoeff2 locCoeff3 locCoeff4 scaleCoeff1
## locCoeff1 9.954e+00 -4.245e-03 1.292e-02 -4.676e-05 1.695e+00
## locCoeff2 -4.245e-03 9.010e-05 -6.018e-05 1.632e-07 -6.348e-04
## locCoeff3 1.292e-02 -6.018e-05 2.533e-04 7.635e-07 2.467e-03
## locCoeff4 -4.676e-05 1.632e-07 7.635e-07 1.660e-08 -2.541e-05
## scaleCoeff1 1.695e+00 -6.348e-04 2.467e-03 -2.541e-05 1.077e+00
## shapeCoeff1 -3.242e-02 8.809e-05 -6.516e-05 2.652e-07 9.445e-04
## shapeCoeff1
## locCoeff1 -3.242e-02
## locCoeff2 8.809e-05
## locCoeff3 -6.516e-05
## locCoeff4 2.652e-07
## scaleCoeff1 9.445e-04
## shapeCoeff1 3.281e-04
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 547“Prediction”
library(fields)
x <- seq(min(coord[,1]), max(coord[,1]), length = 100)
y <- seq(min(coord[,2]), max(coord[,2]), length = 100)
grid <- expand.grid(x, y); colnames(grid) <- c("lon", "lat")
## Watch out we scale the data so we need to use the same transformation
grid[,1] <- grid[, 1] - mean(coord[, 1])
grid[,2] <- grid[, 2] - mean(coord[, 2])
ans <- predict(M0, newdata = grid, ret.per = 20)$Q20
## Do the plot
maps::map(xlim = range(x) / 111, ylim = range(y) / 111)
image(x / 111, y / 111, matrix(ans, 100), add = TRUE, col = tim.colors(64))
contour(x / 111, y / 111, matrix(ans, 100), add = TRUE); maps::map(add = TRUE)
Bayesian Hierarchical Model
Plot the spatial domain and station locations
data(rainfall)
data(swissalt)
par(mar = rep(0, 4), ps = 16)
image(lon.vec, lat.vec, alt.mat, col = terrain.colors(64), asp = 1,
bty = "n", xlab = "n", ylab = "n", axes = FALSE)
swiss(add = TRUE, city = TRUE)
points(coord, pch = 16, cex = 0.75, col = "blue")
Prior specification
hyper <- list()
hyper$betaMeans <- list(loc = rep(0, 3), scale = rep(0, 3), shape = 0)
hyper$betaIcov <- list(loc = diag(rep(1/1000, 3)), scale = diag(rep(1 / 1000, 3)), shape = 1 / 10)
hyper$sills <- list(loc = c(1, 12), scale = c(1, 1), shape = c(1, 0.04))
hyper$ranges <- list(loc = c(5, 3), scale = c(5, 3), shape = c(5, 3))
hyper$smooths <- list(loc = c(1, 1), scale = c(1, 1), shape = c(1, 1))
prop <- list(gev = c(3, 0.1, 0.3), ranges = c(1, 0.8, 1.2), smooths = rep(0, 3))
start <- list(sills = c(10, 10, 0.5), ranges = c(20, 10, 10), smooths = c(1, 1, 1), beta = list(loc = c(25, 0, 0), scale = c(33, 0, 0), shape = 0.001))Running the Gibbs sampler
loc.form <- scale.form <- y ~ lon + lat; shape.form <- y ~ 1
chain <- latent(rain, coord[, 1:2], "powexp", loc.form, scale.form, shape.form,
hyper = hyper, prop = prop,
start = start, n = 100, burn.in = 500, thin = 5)
chain## Effective length: 100
## Burn-in: 500
## Thinning: 5
## Effective NoP: 75.52989
## DIC: 29166.1
##
## Regression Parameters:
## Location Parameters:
## lm1 lm2 lm3
## ci.lower 1.09608 -0.01467 -0.19594
## post.mean 30.50024 0.03854 -0.12421
## ci.upper 64.82665 0.08884 -0.05800
##
## Scale Parameters:
## lm1 lm2 lm3
## ci.lower -5.30524 0.00193 -0.05280
## post.mean 5.95500 0.01644 -0.03197
## ci.upper 17.89916 0.03166 -0.01415
##
## Shape Parameters:
## lm1
## ci.lower 0.06993
## post.mean 0.16673
## ci.upper 0.26259
##
##
## Latent Parameters:
## Location Parameters:
## Covariance family: powexp
## sill range smooth
## ci.lower 5.718 13.520 1.000
## post.mean 9.930 23.864 1.000
## ci.upper 17.044 45.431 1.000
##
## Scale Parameters:
## Covariance family: powexp
## sill range smooth
## ci.lower 0.1274 5.9886 1.0000
## post.mean 0.3122 17.6629 1.0000
## ci.upper 0.5957 27.9217 1.0000
##
## Shape Parameters:
## Covariance family: powexp
## sill range smooth
## ci.lower 0.004151 13.704008 1.000000
## post.mean 0.007827 23.628722 1.000000
## ci.upper 0.013391 36.571825 1.000000Making 20-year return level map
map.latent(chain, ret.per = 20, plot.contour = FALSE, col = tim.colors())
Max-stable processes
USHCN data
data(USHCNTemp)
maps::map("usa")
points(metadata[, c("lon","lat")], pch = 16, cex = 0.75)
F– madogram
data <- maxima.summer
coord <- as.matrix(metadata[, c("lon", "lat")])
fmadogram(data, coord)
fmadogram(data, coord, which = 'ext', n.bins = 100, pch = 16)
abline(h = 1, lty = 2); abline(h = 2, lty = 2)
Simulation from various max-stable models
x <- seq(0, 10, length = 100)
n.sim <- 5
schlather <- rmaxstab(n.sim, x, "powexp", nugget = 0, range = 3, smooth = 1)
extremalt <- rmaxstab(n.sim, x, "tpowexp", DoF = 4, nugget = 0, range = 3, smooth = 1)
brown <- rmaxstab(n.sim, x, "brown", range = 3, smooth = 1)
smith <- rmaxstab(n.sim, x, "gauss", var = 2)
plot(x, schlather[1,], type = "l", ylim = range(schlather))
for (i in 2:5) lines(x, schlather[i,], col = i)
plot(x, smith[1,], type = "l", ylim = range(smith))
for (i in 2:5) lines(x, smith[i,], col = i)
Two-step max-stable process fitting
data(rainfall)
alt <- coord[,3]; coord <- coord[,-3]
## Transformation to unit Frechet margins
rain.frech <- apply(rain, 2, gev2frech, emp = TRUE)
(fit <- fitmaxstab(rain.frech, coord, "powexp", nugget = 0))## Estimator: MPLE
## Model: Schlather
## Weighted: FALSE
## Pair. Deviance: 1136875
## TIC: 1137468
## Covariance Family: Powered Exponential
##
## Estimates
## Marginal Parameters:
## Assuming unit Frechet.
##
## Dependence Parameters:
## range smooth
## 38.4402 0.8528
##
## Standard Errors
## range smooth
## 8.713 0.119
##
## Asymptotic Variance Covariance
## range smooth
## range 75.90902 -0.91383
## smooth -0.91383 0.01417
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 71One-step model fitting
loc.form <- y ~ lon + lat + alt
scale.form <- y ~ lon + lat
shape.form <- y ~ 1
(fit <- fitmaxstab(rain, coord, "powexp", nugget=0, loc.form, scale.form, shape.form,
marg.cov = cbind(alt = alt)))## Computing appropriate starting values
## Starting values are defined
## Starting values are:
## range smooth locCoeff1 locCoeff2 locCoeff3 locCoeff4
## 27.257344723 0.974059614 29.764826111 0.037985247 -0.135302350 0.008372902
## scaleCoeff1 scaleCoeff2 scaleCoeff3 shapeCoeff1
## 10.529521199 0.013058818 -0.040111724 0.156381541## Estimator: MPLE
## Model: Schlather
## Weighted: FALSE
## Pair. Deviance: 2232232
## TIC: 2243533
## Covariance Family: Powered Exponential
##
## Estimates
## Marginal Parameters:
## Location Parameters:
## locCoeff1 locCoeff2 locCoeff3 locCoeff4
## 28.660445 0.039237 -0.134847 0.008425
## Scale Parameters:
## scaleCoeff1 scaleCoeff2 scaleCoeff3
## 8.59677 0.01565 -0.03937
## Shape Parameters:
## shapeCoeff1
## 0.1911
## Dependence Parameters:
## range smooth
## 29.7692 0.9751
##
## Standard Errors
## range smooth locCoeff1 locCoeff2 locCoeff3 locCoeff4
## 9.6080167 0.1744583 6.9550328 0.0097462 0.0166063 0.0005552
## scaleCoeff1 scaleCoeff2 scaleCoeff3 shapeCoeff1
## 4.8427858 0.0071586 0.0115508 0.0518712
##
## Asymptotic Variance Covariance
## range smooth locCoeff1 locCoeff2 locCoeff3
## range 9.231e+01 -1.482e+00 6.902e+00 1.389e-03 -1.886e-02
## smooth -1.482e+00 3.044e-02 -6.333e-02 1.543e-05 9.781e-05
## locCoeff1 6.902e+00 -6.333e-02 4.837e+01 -5.307e-02 -4.118e-02
## locCoeff2 1.389e-03 1.543e-05 -5.307e-02 9.499e-05 -4.612e-05
## locCoeff3 -1.886e-02 9.781e-05 -4.118e-02 -4.612e-05 2.758e-04
## locCoeff4 7.735e-04 -1.761e-05 5.940e-04 -1.975e-06 2.533e-06
## scaleCoeff1 3.595e+00 -5.025e-02 1.111e+01 -1.342e-02 -3.458e-03
## scaleCoeff2 9.077e-03 -5.956e-05 -1.033e-02 2.306e-05 -2.246e-05
## scaleCoeff3 -1.761e-02 8.477e-05 -1.183e-02 -1.005e-05 6.792e-05
## shapeCoeff1 2.473e-01 -2.315e-03 4.184e-02 -2.615e-05 -9.515e-05
## locCoeff4 scaleCoeff1 scaleCoeff2 scaleCoeff3 shapeCoeff1
## range 7.735e-04 3.595e+00 9.077e-03 -1.761e-02 2.473e-01
## smooth -1.761e-05 -5.025e-02 -5.956e-05 8.477e-05 -2.315e-03
## locCoeff1 5.940e-04 1.111e+01 -1.033e-02 -1.183e-02 4.184e-02
## locCoeff2 -1.975e-06 -1.342e-02 2.306e-05 -1.005e-05 -2.615e-05
## locCoeff3 2.533e-06 -3.458e-03 -2.246e-05 6.792e-05 -9.515e-05
## locCoeff4 3.082e-07 -4.400e-04 1.178e-06 -1.332e-06 -6.617e-06
## scaleCoeff1 -4.400e-04 2.345e+01 -2.758e-02 -1.355e-02 4.509e-02
## scaleCoeff2 1.178e-06 -2.758e-02 5.125e-05 -3.135e-05 2.846e-07
## scaleCoeff3 -1.332e-06 -1.355e-02 -3.135e-05 1.334e-04 -1.169e-04
## shapeCoeff1 -6.617e-06 4.509e-02 2.846e-07 -1.169e-04 2.691e-03
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 911plot(fit)
LRT Model selection
M0 <- fitmaxstab(rain, coord, "powexp", nugget = 0, locCoeff4 = 0, loc.form, scale.form,
shape.form, marg.cov = cbind(alt = alt))## Computing appropriate starting values
## Starting values are defined
## Starting values are:
## range smooth locCoeff1 locCoeff2 locCoeff3 locCoeff4
## 27.257344723 0.974059614 29.764826111 0.037985247 -0.135302350 0.008372902
## scaleCoeff1 scaleCoeff2 scaleCoeff3 shapeCoeff1
## 10.529521199 0.013058818 -0.040111724 0.156381541M1 <- fitmaxstab(rain, coord, "powexp", nugget = 0,
loc.form, scale.form, shape.form, marg.cov = cbind(alt = alt))## Computing appropriate starting values
## Starting values are defined
## Starting values are:
## range smooth locCoeff1 locCoeff2 locCoeff3 locCoeff4
## 27.257344723 0.974059614 29.764826111 0.037985247 -0.135302350 0.008372902
## scaleCoeff1 scaleCoeff2 scaleCoeff3 shapeCoeff1
## 10.529521199 0.013058818 -0.040111724 0.156381541anova(M0, M1)## Eigenvalue(s): 161.12
##
## Analysis of Variance Table
## MDf Deviance Df Chisq Pr(> sum lambda Chisq)
## M0 9 2246053
## M1 10 2232232 1 13821 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model selection via TIC
M1 <- fitmaxstab(rain, coord, "powexp", nugget = 0, loc.form, scale.form, shape.form,
marg.cov = cbind(alt = alt))## Computing appropriate starting values
## Starting values are defined
## Starting values are:
## range smooth locCoeff1 locCoeff2 locCoeff3 locCoeff4
## 27.257344723 0.974059614 29.764826111 0.037985247 -0.135302350 0.008372902
## scaleCoeff1 scaleCoeff2 scaleCoeff3 shapeCoeff1
## 10.529521199 0.013058818 -0.040111724 0.156381541M2 <- fitmaxstab(rain, coord, "whitmat", loc.form, scale.form,
shape.form, marg.cov = cbind(alt = alt))## Computing appropriate starting values
## Starting values are defined
## Starting values are:
## nugget range smooth locCoeff1 locCoeff2
## 9.754459e-05 3.511644e+01 3.877117e-01 2.976483e+01 3.798525e-02
## locCoeff3 locCoeff4 scaleCoeff1 scaleCoeff2 scaleCoeff3
## -1.353023e-01 8.372902e-03 1.052952e+01 1.305882e-02 -4.011172e-02
## shapeCoeff1
## 1.563815e-01TIC(M1, M2)## M1 M2
## 2243533 2243665Parameters and return level maps
xlim <- range(M2$coord[,1]); ylim <- range(M2$coord[,2])
data(swissalt)##Retrieve altitude at other places
idx.x <- which(lon.vec >= xlim[1] & lon.vec <= xlim[2])
idx.y <- which(lat.vec >= ylim[1] & lat.vec <= ylim[2])
x <- lon.vec[idx.x]; y <- lat.vec[idx.y]
n.x <- length(x); n.y <- length(y)
covariates <- array(alt.mat[idx.x, idx.y], dim = c(n.x, n.y, 1), dimnames = list(NULL, NULL, "alt"))
covariates[is.na(covariates)] <- 0
par(mar = rep(0, 4), mfrow = c(1, 3))
SpatialExtremes::map(M2, x, y, covariates, param = "loc",
xaxt = "n", yaxt = "n", col = tim.colors()); swiss(add = TRUE)
SpatialExtremes::map(M2, x, y, covariates, param = "scale",
xaxt = "n", yaxt = "n", col = tim.colors()); swiss(add = TRUE)
SpatialExtremes::map(M2, x, y, covariates, param = "quant", ret.per = 20,
xaxt = "n", yaxt = "n", col = tim.colors()); swiss(add = TRUE)